

x <- seq(0,5,0.1)  #横坐标等差数列
y <- df(x,2,19)  #密度函数的值
plot(x,y,type='l') 
x1 <- x[x>=3.7113]
y1 <- y[x>=3.7113]
polygon(c(x1,tail(x1,1),head(x1,1)),c(y1,0,0),col='blue')
abline(h=0)
abline(v=0)

mydata01=juul[c('igf1','tanner')]
cc=complete.cases(mydata01)
mydata02=mydata01[cc,]
nrow(mydata02)

attach(red.cell.folate)
red.cell.folate$rof <- rank(folate)
rofbar <- tapply(rof,ventilation,mean)
ni <- tapply(rof,ventilation,length)
kw <- sum(ni*(rofbar-mean(1:22))^2)*12/22/23


attach(heart.rate)
head(heart.rate)
str(hr)
str(subj)
str(time)
?heart.rate
anova(lm(hr~subj+time))

my.hr <- c(96,110,89,95,128,100,72,79,100,92,106,
            + 86,78,124,98,68,75,106,86,108,85,78,118,100,
            + 67,74,104,92,114,83,83,118,94,71,74,102)
my.subj <- gl(9,1,36)
my.time <- gl(4,9,36,labels=c(0,30,60,90))
my.heart.rate <- data.frame(my.hr, my.subj, my.time)
my.heart.rate


attach(thuesen)
lm01 <- lm(short.velocity~blood.glucose)
summary(lm01)
anova(lm01)

walk <- unlist(zelazo)
group <- factor(rep(1:4,c(6,6,6,5)), labels=names(zelazo))
summary(lm(walk ~ group))
t.test(zelazo$active,zelazo$ctr.8w) # first vs. last
t.test(zelazo$active,unlist(zelazo[-1])) # first vs. rest

library(ISwR)  #1 
zelazo  #2 
?zelazo  #3 
walk <- unlist(zelazo)  #4 
group <- factor(rep(1:4,c(6,6,6,5)),labels=names(zelazo))  #5 
mydata <- data.frame(walk=walk,group=group)  #6 
mydata

x1 <- zelazo$active
x2 <- zelazo$passive
x3 <- zelazo$none
x4 <- zelazo$ctr.8w
x <- c(x1,x2,x3,x4)
x1bar <- mean(x1)
x2bar <- mean(x2)
x3bar <- mean(x3)
x4bar <- mean(x4)
xbar <- mean(x)
n1 <- length(x1)
n2 <- length(x2)
n3 <- length(x3)
n4 <- length(x4)

SSDw <- sum((x1-x1bar)^2)+sum((x2-x2bar)^2)+sum((x3-x3bar)^2)+sum((x4-x4bar)^2)
SSDb <- n1*(x1bar-xbar)^2+n2*(x2bar-xbar)^2+n3*(x3bar-xbar)^2+n4*(x4bar-xbar)^2
SSDw
SSDb

N=length(x)
k=4
MSw <- SSDw/(N-k)
MSb <- SSDb/(k-1)
myf <- MSb/MSw
1-pf(myf,3,19)

lm01 <- lm(walk~group,data=mydata)
summary(lm01)
anova(lm01)

pairwise.t.test(walk,group)
pairwise.t.test(walk,group,p.adj='bonferroni')
t.test(zelazo$active,zelazo$ctr.8w)
oneway.test(walk~group)
bartlett.test(walk~group)
kruskal.test(walk~group)

walkbar <- tapply(walk,group,mean)  #1 
walksd <- tapply(walk,group,sd)  #2 
walkn <- tapply(walk,group,length)  #3 
sem <- walksd/sqrt(walkn)  #4 
stripchart(walk~group,method='jitter',jitter=0.05,pch=16,vert=T)  #5 
arrows(1:4,walkbar+sem,1:4,walkbar-sem,angle=45,code=3,length=0.1)  #6 
lines(1:4,walkbar,pch=1,type='b',cex=1)  #7 





lung
?lung
lm07 <- lm(volume~method+subject,data=lung)
summary(lm07)


lung
?lung
fit <- lm(volume~method+subject, data=lung)
anova(fit)
summary(fit)


attach(lung)
interaction.plot(method,subject,volume)

attach(lung)
x <- lung$volume
xbar <- mean(x)
xidotbar <- tapply(volume,subject,mean) #每行的平均值
xdotjbar <- tapply(volume,method,mean) #每列的平均值
m <- 6 #行数
n <- 3 #列数
SSDr <- n*sum((xidotbar-xbar)^2) #行间方差
SSDc <- m*sum((xdotjbar-xbar)^2) #列间方差
SSDt <- sum((x-xbar)^2) #总方差
SSDres <- SSDt-SSDr-SSDc #按公式计算残差方差
xi <- rep(xidotbar,each=3) 
xj <- rep(xdotjbar,6)
mydata <- lung
mydata$xi <- xi
mydata$xj <- xj
SSDres02 <- sum((x-xi-xj+xbar)^2) #按定义计算残差方差

attach(lung)
lm01 <- lm(volume ~ method + subject)
summary(lm01)
anova(lm01)
friedman.test(volume ~ subject | method)
friedman.test(volume ~ method | subject)

myfr <- SSDr/(m-1)/SSDres*(m-1)*(n-1)
1-pf(myfr,m-1,(m-1)*(n-1))
myfc <- SSDc/(n-1)/SSDres*(m-1)*(n-1)
1-pf(myfc,n-1,(m-1)*(n-1))


walk <- unlist(zelazo)
group <- factor(rep(1:4,c(6,6,6,5)), labels=names(zelazo))

kruskal.test(walk ~ group)
wilcox.test(zelazo$active,zelazo$ctr.8w) # first vs. last
wilcox.test(zelazo$active,unlist(zelazo[-1])) # first vs. rest
friedman.test(volume ~ method | subject, data=lung)
wilcox.test(lung$volume[lung$method=="A"],
            lung$volume[lung$method=="C"], paired=TRUE) # etc.


attach(juul)
tapply(sqrt(igf1),tanner, sd, na.rm=TRUE)
plot(sqrt(igf1)~jitter(tanner))
oneway.test(sqrt(igf1)~tanner)
pairwise.t.test(sqrt(igf1),tanner, pool.sd = F)




tanner <- factor(tanner,labels=c('i','ii','iii','iv','v'))
lm01 <- lm(igf1~tanner)
summary(lm01)
anova(lm01)




library(pracma)
A=matrix(c(3,-2,0,-1,3,-1,-5,7,-1),nrow=3,byrow=T);print(A);eig(A)
B=matrix(c(4,-5,7,1,-4,9,-4,0,5),nrow=3,byrow=T);print(B);eig(B)
C=matrix(c(3,6,6,0,2,0,-3,-12,-6),nrow=3,byrow=T);print(C);eig(C)

A=matrix(c(3,4,2,0,2,3,2,3),nrow=2,byrow=T);print(A);rref(A)

A=matrix(c(2,0,2,3),nrow=2,byrow=T)
P=matrix(c(3,4,2,3),nrow=2,byrow=T)
print(inv(P)%*%A%*%P)
X=matrix(c(1,2),nrow=2,byrow=T)
print(inv(P)%*%A%*%X)


x1=c(1073, 1009, 1060, 1001, 1002, 1012, 1009, 1028)
x2=c(1107, 1092,  990, 1109, 1090, 1074, 1122, 1001)
x3=c(1093, 1029, 1080, 1021, 1022, 1032, 1029, 1048)
x=c(x1,x2,x3)
A=gl(3,8,24)
levels(A)=c('A1','A2','A3')
mydata=data.frame(x,A)
mydata
anova(lm(x~A,data=mydata))



x <- seq(0,5,0.05)  #横坐标等差数列
y <- df(x,5,12)  #密度函数的值
plot(x,y,type='l',ylab='y = df(x,5,12)',ylim=c(0,0.8))  #画出密度函数的图像
x1 <- x[x>=2.5]  # F统计值及其右边（阴影部分的下界）
y1 <- y[x>=2.5]  # 阴影部分的上界
polygon(c(x1,tail(x1,1),head(x1,1)),c(y1,0,0),col='blue')
abline(h=0)  #横坐标轴
abline(v=0)  #纵坐标轴




